function PDE = process_pde_t_2d(example_index)
switch example_index
    case 1
        nu = 1;
        PDE.index = example_index;
        PDE.nu = nu;
        PDE.u1 = @(x,y,t) (x.^2.*y.^2 + exp(-y)).*cos(2*pi*t);
        PDE.u2 = @(x,y,t) (-(2/3)*x.*y.^3 + 2 - pi*sin(pi*x)).*cos(2*pi*t);
        PDE.p = @(x,y,t) -(2-pi*sin(pi*x)).*cos(2*pi*y).*cos(2*pi*t);
        PDE.u1dx = @(x,y,t) 2*x.*y.^2.*cos(2*pi*t);
        PDE.u1dy = @(x,y,t) -cos(2*pi*t).*(-2*y.*x.^2 + exp(-y));
        PDE.u2dx = @(x,y,t) -cos(2*pi*t).*((2*y.^3)./3 + pi^2*cos(x*pi));
        PDE.u2dy = @(x,y,t) -2*x.*y.^2.*cos(2*pi*t);
        PDE.pdx = @(x,y,t) pi^2*cos(2*pi*t).*cos(pi*x).*cos(2*pi*y);
        PDE.pdy = @(x,y,t) -2*pi*cos(2*pi*t).*sin(2*pi*y).*(pi*sin(pi*x) - 2);
        PDE.f1 = @(x,y,t) -2*pi*(x.^2.*y.^2 + exp(-y)).*sin(2*pi*t) + (-2*nu*x.^2 - 2*nu*y.^2 - nu*exp(-y) + pi^2*cos(pi*x).*cos(2*pi*y)).*cos(2*pi*t);
        PDE.f2 = @(x,y,t) -2*pi*(-(2/3)*x.*y.^3 + 2 - pi*sin(pi*x)).*sin(2*pi*t) + (4*nu*x.*y - nu*pi^3*sin(pi*x) + 2*pi*(2 - pi*sin(pi*x)).*sin(2*pi*y)).*cos(2*pi*t);
        PDE.domain = [0, 1, -0.25, 0, 0, 1];
        PDE.bdary = [0, 0, 0, 0];
    case 2
        PDE.index = example_index;
        nu = 1;
        PDE.nu = nu;
        PDE.u1 = @(x,y,t) pi*sin(2*pi*y).*sin(pi*x).^2.*sin(t);
        PDE.u2 = @(x,y,t) -pi*sin(2*pi*x).*sin(pi*y).^2.*sin(t);
        PDE.p = @(x,y,t) 10*cos(pi*x).*sin(pi*y).*sin(t);
        PDE.u1dx = @(x,y,t) 2*pi^2*cos(pi*x).*sin(pi*x).*sin(2*pi*y).*sin(t);
        PDE.u1dy = @(x,y,t) 2*pi^2*cos(2*pi*y).*sin(pi*x).^2.*sin(t);
        PDE.u2dx = @(x,y,t) -2*pi^2*cos(2*pi*x).*sin(pi*y).^2.*sin(t);
        PDE.u2dy = @(x,y,t) -2*pi^2*cos(pi*y).*sin(2*pi*x).*sin(pi*y).*sin(t);
        PDE.pdx = @(x,y,t) -10*pi*sin(pi*x).*sin(pi*y).*sin(t);
        PDE.pdy = @(x,y,t) 10*pi*cos(pi*x).*cos(pi*y).*sin(t);
        PDE.f1 = @(x,y,t) pi*sin(pi*x).^2.*sin(2*pi*y).*cos(t) - nu*(2*pi^3*cos(pi*x).^2.*sin(2*pi*y).*sin(t) - 6*pi^3*sin(pi*x).^2.*sin(2*pi*y).*sin(t)) - 10*pi*sin(pi*x).*sin(pi*y).*sin(t);
        PDE.f2 = @(x,y,t) nu*(2*pi^3*cos(pi*y).^2.*sin(2*pi*x).*sin(t) - 6*pi^3*sin(2*pi*x).*sin(pi*y).^2.*sin(t)) - pi*sin(2*pi*x).*sin(pi*y).^2.*cos(t) + 10*pi*cos(pi*x).*cos(pi*y).*sin(t);
        PDE.domain = [0, 1, 0, 1, 0, 1];
        PDE.bdary = [0, 0, 0, 0];
    case 3
        PDE.index = example_index;
        nu = 1;
        PDE.nu = nu;
        PDE.u1 = @(x,y,t) 10.*x.^2.*y.*cos(t).*(2.*y - 1).*(x - 1).^2.*(y - 1);
        PDE.u2 = @(x,y,t) -10.*x.*y.^2.*cos(t).*(2.*x - 1).*(x - 1).*(y - 1).^2;
        PDE.p = @(x,y,t) cos(t).*(20.*x - 10).*(2.*y - 1);
        PDE.u1dx = @(x,y,t) 20*x.*y.*cos(t).*(2*y - 1).*(x - 1).^2.*(y - 1) + 10*x.^2.*y.*cos(t).*(2*x - 2).*(2*y - 1).*(y - 1);
        PDE.u1dy = @(x,y,t) 10*x.^2.*cos(t).*(2*y - 1).*(x - 1).^2.*(y - 1) + 20*x.^2.*y.*cos(t).*(x - 1).^2.*(y - 1) + 10*x.^2.*y.*cos(t).*(2*y - 1).*(x - 1).^2;
        PDE.u2dx = @(x,y,t) -10*y.^2.*cos(t).*(2*x - 1).*(x - 1).*(y - 1).^2 - 20*x.*y.^2.*cos(t).*(x - 1).*(y - 1).^2 - 10*x.*y.^2.*cos(t).*(2*x - 1).*(y - 1).^2;
        PDE.u2dy = @(x,y,t) -20.*x.*y.*cos(t).*(2*x - 1).*(x - 1).*(y - 1).^2 - 10*x.*y.^2.*cos(t).*(2*x - 1).*(2*y - 2).*(x - 1);
        PDE.pdx = @(x,y,t) 20.*cos(t).*(2*y - 1);
        PDE.pdy = @(x,y,t) 2*cos(t).*(20*x - 10);
        PDE.f1 = @(x,y,t) 20*cos(t).*(2*y - 1) - nu*(40*x.^2.*cos(t).*(x - 1).^2.*(y - 1) + 20*x.^2.*cos(t).*(2*y - 1).*(x - 1).^2 + 40*x.^2.*y.*cos(t).*(x - 1).^2 + 20*x.^2.*y.*cos(t).*(2*y - 1).*(y - 1) + 20*y.*cos(t).*(2*y - 1).*(x - 1).^2.*(y - 1) + 40*x.*y.*cos(t).*(2*x - 2).*(2*y - 1).*(y - 1)) - 10*x.^2.*y.*sin(t).*(2*y - 1).*(x - 1).^2.*(y - 1);
        PDE.f2 = @(x,y,t) 2*cos(t).*(20*x - 10) + nu*(40*y.^2.*cos(t).*(x - 1).*(y - 1).^2 + 20*y.^2.*cos(t).*(2*x - 1).*(y - 1).^2 + 40*x.*y.^2.*cos(t).*(y - 1).^2 + 20*x.*y.^2.*cos(t).*(2*x - 1).*(x - 1) + 20*x.*cos(t).*(2*x - 1).*(x - 1).*(y - 1).^2 + 40*x.*y.*cos(t).*(2*x - 1).*(2*y - 2).*(x - 1)) + 10*x.*y.^2.*sin(t).*(2*x - 1).*(x - 1).*(y - 1).^2;
        PDE.domain = [0, 1, 0, 1, 0, 1];
        PDE.bdary = [0, 0, 0, 0];
    otherwise
        error("Invalid PDE index.");
end
fprintf("PDE index: %d\n", PDE.index);
fprintf("PDE domain: xmin,xmax,ymin,ymax,tmin,tmax\n");
disp(PDE.domain);
fprintf("PDE boundary: bottom,right,top,left (0=Dirichlet,1=Neumann,2=Robin)\n");
disp(PDE.bdary);
end